Librairies

rm(list=objects())
graphics.off()

library(INLA)
library(inlabru)
library(sp)
library(spatstat)
library(ggplot2)
library(raster)       # raster altitude
library(numbers)      # fonction mod
library(viridis)      # couleurs graphique altitude
library(gstat)        # krigeage
library(dbmss)        # fonction Mhat

Chargement des données

Les données fournies sont celles de la parcelle 16 de Paracou (année 2015) auxquelles j’ai rajouté 3 variables.

• X et Y : les coordonnées sur la parcelle, avec pour origine le coin Sud-Ouest et pour axes les côtés Sud et Ouest du carré,
• altitude : l’altitude du terrain à l’emplacement de l’arbre.

J’ai également supprimé les arbres qui sortaient de la parcelle. Pour plus de précisions concernant la construction de ces données, voir le fichier “preparation_donnees_paracou.R”.

load("data/p16_avec coord_et_altitude_2015.RData")

Fonction de préparation des données

Cette fonction crée à partir des coordonnées x et y du semis de points, tous les objets dont nous aurons besoin pour notre étude. Elle retourne le nombre d’observations (nb_obs), un SpatialPointsDataFrame de coordonnées (coord_pts), un objet ppp des coordonnées (coord_pts_ppp), un DataFrame des coordonnées (df) et un mesh (mesh) adapté aux coordonnées.

L’argument data doit être composé de deux colonnes nommées x et y. L’argument mesh = TRUE par défaut. L’argument size spécifie la taille d’un côté du domaine (par défaut = 5), et donc les valeurs max. des coordonnées entrées en argument (deux valeurs utilisées : 5 ou 500).

Sûrement du à un problème d’implémentation, nous ne sommes pas en mesure de construire un maillage régulier sur des données allant de 0 à 500 (valeurs trop élevées qui engendre un temps de calcul démesurément long). L’argument size fixé à une valeur de 500 servira donc uniquement à tracer certains graphiques en taille réelle ; si size = 500, la fonction ne produira jamais de maillage. La solution provisoire trouvée pour nos études est de considérer un domaine de taille 5x5, et de diviser les coordonnées des arbres de la parcelle par 100.

preparation = function(data, mesh = TRUE, size = 5){
  
sortie = list()
  
# spatialPointDataFrame  
coord_pts = data.frame(data)                   
sortie$nb_obs = nrow(coord_pts)
coordinates(coord_pts) = c("x", "y")           
sortie$coord_pts = coord_pts

# dataframe pour les graphiques 
sortie$df = data.frame(data)

# création de d'objet ppp point pattern à partir des coordonnées
sortie$coord_pts_ppp = ppp(coord_pts$x, coord_pts$y, c(0,size), c(0,size))


if(size==5)
{  if(mesh==TRUE)
  {   # mesh apadté aux données
      bnd = spoly(data.frame(lon =  c(0, 5, 5, 0, 0), lat = c(0, 0, 5, 5, 0))) 
      boundary <- list(as.inla.mesh.segment(bnd), NULL)
      mesh.loc <- coord_pts
      mesh <- inla.mesh.2d(loc=mesh.loc,
                     boundary=boundary,
                     max.edge=c(0.2, 0.74),
                     min.angle=c(30, 21),
                     max.n=c(48000, 16000),
                     max.n.strict=c(128000, 128000), 
                     cutoff=0.008,
                     offset=c(0.16, 0.43))
      sortie$mesh = mesh
  }
}

return(sortie)}

Mesh régulier, frontières et graph

Comme expliqué précedemment, les maillages réguliers sur des domaines dont les côtés dépassent 50 sont extrêmement longs à construire, même avec des exigences de qualité très faibles… La solution provisoire que nous allons donc adopter pour ces études est d’effectuer notre analyse sur un domaine de taille 5X5 et de diviser les coordonnées de chaque arbre par 100.

# frontière du domaine 5x5 (sens inverse des aiguilles d'une montre)
bnd = spoly(data.frame(lon =  c(0, 5, 5, 0, 0), lat = c(0, 0, 5, 5, 0)))  

# même frontière en objet owin
bnd_owin = owin(xrange=c(0,5), yrange=c(0,5))

# mesh régulier sur le domaine 5x5
mesh_r = inla.mesh.2d(boundary=bnd, max.edge = 0.2) 
ggplot() + gg(mesh_r) + coord_fixed() + ggtitle("Mesh régulier")

# grille utilisée pour les prédictions sur le domaine
pix = pixels(mesh_r, nx = 200, ny = 200)   

# coordonnées des noeuds de cette grille dans un dataframe
pix_df = as.data.frame(pix)


# fonction qui définit le min et le max de l'échelle de couleurs pour l'affichage
colsc <- function(max, min=0) {scale_fill_gradientn(colours = rev(RColorBrewer::brewer.pal(11,"RdYlBu")), limits = range(c(min, max), na.rm=TRUE))}

Poisson homogène

Pour effectuer une régression de Poisson homogène, la fonction bru se base sur une grille arbitraire couvrant l’ensemble du domaine. Elle prend comme argument un SpatialPointsDataFrame contenant (pour coordonnées) les centres des cellules de cette grille et (pour valeur) le nombre de points du semis contenus dans chaque cellule correspondante. Il nous faut donc au préalable construire cette grille et son DataFrame.

# Simulation
poish = preparation(rpoispp(lambda = 10, win=bnd_owin), mesh = FALSE)


# Grille pour la régression de poisson
nb_carres_par_cote = 10                                                # résolution de la grille (arbitraire)
grid = raster(xmn=0, ymn=0, xmx=5, ymx=5, res=5/nb_carres_par_cote)    # création de la grille sur le domaine
grid[] = 0                                                             # grid[] : nombre d'observations dans chaque cellule

tab = table(cellFromXY(grid, poish$coord_pts))          # compte pour chaque cellule de la grille, le nombre d'observations  
grid[as.numeric(names(tab))] = tab                      # qu'elle contient,
count_df = data.frame(coordinates(grid), count=grid[])  # ce résultat est ensuite stocké dans la variable count 
coordinates(count_df) <- c("x", "y")                    # conversion en SpatialPointsDataFrame
  

# Régression
fit_poish = bru(count ~ Intercept, count_df, family="poisson")

lambda_poish = predict(fit_poish, pix, ~ exp(Intercept))
lambda_poish_df = as.data.frame(lambda_poish)                 
# Pour chaque point de pix (en ligne), on a la loi du lambda associé.
lambda_poish_df$intensity = lambda_poish$mean*(nb_carres_par_cote**2/25)          
# Renormalisation pour avoir l'intensité, au lieu du nombre de points par cellule de la grille.

# Plot
ggplot() + gg(poish$coord_pts) + gg(bnd) + coord_fixed()

## Regions defined for each Polygons

plot(grid); points(poish$coord_pts, pch=20)

# Valeurs
print(paste("Intensité théorique 10")) 

## [1] "Intensité théorique 10"

print(paste("Intensité observée", poish$nb_obs/25))             

## [1] "Intensité observée 9.76"

print(paste("Intensité poisson homogène", exp(fit_poish$summary.fixed$mean)*(nb_carres_par_cote**2/25) ))

## [1] "Intensité poisson homogène 9.75994559179058"

summary(fit_poish)

## 
## --- Likelihoods ----------------------------------------------------------------------------------
## 
## 
## --- Criteria -------------------------------------------------------------------------------------
## 
## Watanabe-Akaike information criterion (WAIC):    3.899e+02
## Deviance Information Criterion (DIC):        3.896e+02
## 
## --- Fixed effects -------------------------------------------------------------------------------- 
## 
##                mean         sd 0.025quant  0.5quant 0.975quant      mode
## Intercept 0.8919925 0.06407409  0.7638826 0.8927921   1.015553 0.8943809
##           signif
## Intercept   TRUE
## 
## --- Random effects ------------------------------------------------------------------------------- 
## 
## None.
## 
## --- All hyper parameters (internal representation) ----------------------------------------------- 
## 
## NULL

Poisson non-homogène

Nous souhaitons générer un processus non-homogène sur notre domaine. Poour cela nous choisissons arbitrairement comme intensité théorique sur le champ, la fonction qui somme les coordonées.

# Intensité théorique
intensity_pois = function(x,y){x+y}   
grille = expand.grid(x=seq(0,5,length=200), y=seq(0,5,length=200))
grille$intensity = grille$x + grille$y


# Simulation
pois = preparation(rpoispp(lambda = intensity_pois, win=bnd_owin))

# Intensité observée
densi = density(pois$coord_pts_ppp, sigma=bw.diggle(pois$coord_pts_ppp))


# Echelle de couleurs commune
intensity_max <- ceiling(max(max(densi), max(grille$intensity)))
col_fixe <- colourmap(hcl.colors(intensity_max), breaks=0:intensity_max)

# Plot
ggplot() + gg(pois$coord_pts) + gg(bnd) + coord_fixed()

## Regions defined for each Polygons

mat_pois = matrix(grille$intensity, ncol=200, nrow=200, byrow = TRUE)
plot(im(mat_pois, xcol=seq(ncol(mat_pois)), yrow=seq(nrow(mat_pois))), main="Intensité theorique", col=col_fixe)

plot(densi, main="Intensité observée", col=col_fixe)

Cox homogène

Nous simulons ci-dessous un processus de Cox homogène muni d’une fonction de covariance de Matern, dont les paramètres sont \(moy. = 10\), \(var. = \exp(0.2)\), \(\alpha=1/2\) et \(\nu=1\).

# Simulation
lgcph = rLGCP(model='matern', mu = log(10), var=0.2, scale=1/2, nu=1, win=bnd_owin)
coxh = preparation(data.frame(x = lgcph$y, y = lgcph$x))                 


# lambda_theo est l'intensité exacte de notre processus et c'est une réalisation du champ aléatoire théorique de moyenne 10
lambda_theoh <- attr(lgcph, 'Lambda')
lambda_theoh_df = data.frame(expand.grid(x = lambda_theoh$yrow, y = lambda_theoh$xcol), val=as.vector(lambda_theoh$v))


# Régression
fit_coxh = lgcp(coordinates ~ Intercept, coxh$coord_pts, samplers = bnd)

# Plot
ggplot() + gg(coxh$coord_pts) + gg(bnd) + coord_fixed()

## Regions defined for each Polygons

ggplot(lambda_theoh_df, aes(x=x, y=y)) + geom_raster(aes(fill = val, interpolate=TRUE)) + coord_fixed() + colsc(max(lambda_theoh_df$val)) + gg(coxh$coord_pts) + ggtitle(label = "Lambda théorique")

## Warning: Ignoring unknown aesthetics: interpolate

# Valeurs 
print(paste("Intensité théorique 10"))  

## [1] "Intensité théorique 10"

print(paste("Intensité observée", coxh$nb_obs/25))             

## [1] "Intensité observée 9.64"

print(paste("Intensité cox homogène", exp(fit_coxh$summary.fixed$mean)))

## [1] "Intensité cox homogène 9.70222535866406"

summary(fit_coxh)

## 
## --- Likelihoods ----------------------------------------------------------------------------------
## 
## 
## --- Criteria -------------------------------------------------------------------------------------
## 
## Watanabe-Akaike information criterion (WAIC):    -6.103e+02
## Deviance Information Criterion (DIC):        -6.113e+02
## 
## --- Fixed effects -------------------------------------------------------------------------------- 
## 
##               mean        sd 0.025quant 0.5quant 0.975quant     mode
## Intercept 2.272355 0.0644722   2.143435 2.273165    2.39667 2.274774
##           signif
## Intercept   TRUE
## 
## --- Random effects ------------------------------------------------------------------------------- 
## 
## None.
## 
## --- All hyper parameters (internal representation) ----------------------------------------------- 
## 
## NULL

Cox non-homogène

Sur un processus de Cox non-homogène

Cette fois, nous simulons un processus de Cox non-homogène dont la moyenne n’est alors plus constante (nous choisissons la fonction damier de densité 5 et 20). Le processus est toujours muni d’une fonction de covariance de Matern dont les paramètres sont \(var. = \exp(0.2)\), \(\alpha=1/2\) et \(\nu=1\).

Afin d’effectuer la régression de Cox non-homogène, nous devons définir une variable de type champ aléatoire qui représentera les interactions intra-spécifique et nous permettra d’expliquer l’intensité. Le package INLA apporte une solution inattendue pour cette variable ; ce champ peut en effet être nommé comme on le souhaite et est reconnu seulement grâce à ces arguments (map = …., model = ….). Cette fonction est appelée “f” est dans la documentation R. Dans les exemples ci-dessous elle sera appelée “field”. Nous définissons “field” comme un champ aléatoire de covariance de Matérn, nous définissons donc un sous-modèle avec les paramètres du champ comme inconnus grâce à la fonction inla.spde2.pcmatern.

# Intensité exacte
fct_damier = function(x,y){n = length(x); return(log(5+15*(mod(floor(x)+floor(y), rep(2,n)))))}
grille = expand.grid(x=seq(0,5,length=200), y=seq(0,5,length=200))
grille$intensity = exp(fct_damier(grille$x, grille$y))

# Simulation
lgcp_sim = rLGCP(model='matern', mu=fct_damier, var=0.2, scale=1/2, nu=1, win=bnd_owin)
cox = preparation(data.frame(x = lgcp_sim$y, y = lgcp_sim$x))

# lambda_theo est l'intensité exacte de notre processus et c'est une réalisation du champ aléatoire théorique de moyenne 10
lambda_theo <- attr(lgcp_sim, 'Lambda')
lambda_theo_df = data.frame(expand.grid(x = lambda_theo$yrow, y = lambda_theo$xcol), val=as.vector(lambda_theo$v))

# Régression
model_matern = inla.spde2.pcmatern(cox$mesh, prior.sigma = c(0.1, 0.01), prior.range = c(5, 0.01))
fit_cox = lgcp(coordinates ~ field(map = coordinates, model = model_matern) + Intercept, cox$coord_pts, samplers = bnd)

lambda_cox = predict(fit_cox, pix, ~ exp(field + Intercept))
lambda_cox_df = as.data.frame(lambda_cox)




# Plot
ggplot() + gg(cox$coord_pts) + gg(bnd) + coord_fixed()

## Regions defined for each Polygons

ggplot(grille, aes(x=x, y=y)) + geom_raster(aes(fill = intensity, interpolate=TRUE)) + coord_fixed() +  colsc(max = max(lambda_theo_df$val)) + ggtitle(label = "Intensité théorique")

## Warning: Ignoring unknown aesthetics: interpolate

ggplot(lambda_theo_df, aes(x=x, y=y)) + geom_raster(aes(fill = val, interpolate=TRUE)) + coord_fixed() + colsc(max = max(lambda_theo_df$val)) + gg(cox$coord_pts) + ggtitle(label = "Intensité de la réalisation du champ")

## Warning: Ignoring unknown aesthetics: interpolate

plot(density(cox$coord_pts_ppp, sigma = bw.diggle(cox$coord_pts_ppp)), main="Intensité observée")  

mat_cox = matrix(lambda_cox_df$mean, ncol=200, nrow=200, byrow = TRUE)
plot(im(mat_cox, xcol=seq(ncol(mat_cox)), yrow=seq(nrow(mat_cox))), main="Intensité moyenne") 

mat_cox = matrix(lambda_cox_df$sd, ncol=200, nrow=200, byrow = TRUE)
plot(im(mat_cox, xcol=seq(ncol(mat_cox)), yrow=seq(nrow(mat_cox))), main="Ecart-type") 

summary(fit_cox)

## 
## --- Likelihoods ----------------------------------------------------------------------------------
## 
## 
## --- Criteria -------------------------------------------------------------------------------------
## 
## Watanabe-Akaike information criterion (WAIC):    -8.989e+02
## Deviance Information Criterion (DIC):        -8.991e+02
## 
## --- Fixed effects -------------------------------------------------------------------------------- 
## 
##               mean        sd 0.025quant 0.5quant 0.975quant     mode
## Intercept 2.489473 0.5773392   1.381838 2.489828   3.588412 2.492595
##           signif
## Intercept   TRUE
## 
## --- Random effects ------------------------------------------------------------------------------- 
## 
## field ranges: mean = [-0.00851155, 0.00603428], sd = [0.526497, 0.563169], quantiles = [-1.06471 : 1.06777]
## 
## --- All hyper parameters (internal representation) ----------------------------------------------- 
## 
##                        mean          sd  0.025quant    0.5quant 0.975quant
## Range for field 42.32952059 58.69558380 4.020955550 24.89737502 188.400863
## Stdev for field  0.03129879  0.03642442 0.002680121  0.02035396   0.126565
##                         mode
## Range for field 10.179032859
## Stdev for field  0.007305385
## 
## 
## --- Field 'field' transformed hyper parameters ---
##          param         mean          var           sd            lq
## 1        range 41.848345415 2.871187e+03 53.583455583  3.893303e+00
## 2    log.range  3.244986768 9.562635e-01  0.977887256  1.399332e+00
## 3     variance  0.002147927 4.871517e-05  0.006979626 -6.087022e-05
## 4 log.variance -7.848836246 3.848233e+00  1.961691335 -1.183507e+01
##          median           uq
## 1 24.6982490732 185.14378766
## 2  3.2115729110   5.22760208
## 3  0.0003483903   0.01525113
## 4 -7.7952687180  -4.15498732

Sur un processus de Matérn

La fonction qui simule un processus de Matérn est rMatClust(kappa, scale, mu, win, …. ) de spatstat avec :

kappa = intensité du premier processus de poisson (celui qui génère les centres des cercles) scale = rayon des cercles mu = nombre de points moyen par cercle win = fenêtre globale

# Simulation
matern = rMatClust(kappa=0.3, scale=0.4, mu=20, win = bnd_owin)
mat = preparation(data.frame(x = matern$x, y = matern$y))                 

# Régression
model_matern = inla.spde2.pcmatern(mat$mesh, prior.sigma = c(0.1, 0.01), prior.range = c(5, 0.01))
fit_mat = lgcp(coordinates ~ field(map = coordinates, model = model_matern) + Intercept, mat$coord_pts, samplers = bnd)

lambda_mat = predict(fit_mat, pix, ~ exp(field + Intercept))
lambda_mat_df = as.data.frame(lambda_mat)


# Plot
ggplot() + gg(mat$coord_pts) + gg(bnd) + coord_fixed()

## Regions defined for each Polygons

plot(density(mat$coord_pts_ppp, sigma = bw.diggle(mat$coord_pts_ppp)), main="Intensité observée")  

mat_mat = matrix(lambda_mat_df$mean, ncol=200, nrow=200, byrow = TRUE)
plot(im(mat_mat, xcol=seq(ncol(mat_mat)), yrow=seq(nrow(mat_mat))), main="Intensité moyenne") 

mat_sd = matrix(lambda_mat_df$sd, ncol=200, nrow=200, byrow = TRUE)
plot(im(mat_sd, xcol=seq(ncol(mat_sd)), yrow=seq(nrow(mat_sd))), main="Ecart-type")     

summary(fit_mat)

## 
## --- Likelihoods ----------------------------------------------------------------------------------
## 
## 
## --- Criteria -------------------------------------------------------------------------------------
## 
## Watanabe-Akaike information criterion (WAIC):    -9.579e+02
## Deviance Information Criterion (DIC):        -9.598e+02
## 
## --- Fixed effects -------------------------------------------------------------------------------- 
## 
##               mean        sd 0.025quant 0.5quant 0.975quant     mode
## Intercept 1.439986 0.4833149  0.4907232 1.440095   2.387725 1.440355
##           signif
## Intercept   TRUE
## 
## --- Random effects ------------------------------------------------------------------------------- 
## 
## field ranges: mean = [-2.07012, 3.11968], sd = [0.553607, 1.6963], quantiles = [-5.21758 : 4.46431]
## 
## --- All hyper parameters (internal representation) ----------------------------------------------- 
## 
##                       mean         sd  0.025quant    0.5quant  0.975quant
## Range for field 40.4884204 60.4016291 3.106282508 22.66336008 188.5357188
## Stdev for field  0.0418197  0.0620962 0.002669162  0.02331727   0.1959196
##                        mode
## Range for field 8.041174555
## Stdev for field 0.007020902
## 
## 
## --- Field 'field' transformed hyper parameters ---
##          param         mean          var          sd            lq
## 1        range 39.969590947 2.982100e+03 54.60860853   2.955772819
## 2    log.range  3.140773192 1.087530e+00  1.04284728   1.141436642
## 3     variance  0.004894685 4.806958e-04  0.02192478  -0.007111195
## 4 log.variance -7.528669521 4.758744e+00  2.18145466 -11.837313864
##         median          uq
## 1 22.457889591 185.0392461
## 2  3.117215524   5.2268952
## 3  0.000016481   0.0274249
## 4 -7.523275155  -3.2879265

Représentation de l’altitude

# Altitude pour tous les arbres
altitude = data.frame(x=p16$X,y=p16$Y,z=p16$altitude)
coordinates(altitude)=c("x","y")
ggplot(p16, aes(x=X, y=Y)) + geom_point(aes(color = altitude)) + scale_color_viridis(option="D") + coord_fixed() + ggtitle("Altitude exacte pour tous les arbres")

# Altitude sur un échantillon d'arbres dans la parcelle (nécessaire à la rapidité du krigeage)
diviseur = 10           
sample_p16 = p16[sample(1:nrow(p16), floor(nrow(p16)/diviseur)),]
ggplot(sample_p16, aes(x=X, y=Y)) + geom_point(aes(color = altitude)) + scale_color_viridis(option="D") + coord_fixed() + ggtitle("Altitude exacte sur un échantillon des arbres")

# Visualisation sur la carte
library(rgdal)

## rgdal: version: 1.4-3, (SVN revision 828)
##  Geospatial Data Abstraction Library extensions to R successfully loaded
##  Loaded GDAL runtime: GDAL 2.2.3, released 2017/11/20
##  Path to GDAL shared files: C:/Users/lena.klay/Documents/R/win-library/3.6/rgdal/gdal
##  GDAL binary built with GEOS: TRUE 
##  Loaded PROJ.4 runtime: Rel. 4.9.3, 15 August 2016, [PJ_VERSION: 493]
##  Path to PROJ.4 shared files: C:/Users/lena.klay/Documents/R/win-library/3.6/rgdal/proj
##  Linking to sp version: 1.3-1

library(leaflet)
library(mapview)

coord_utm = data.frame(long=p16$Xutm, lat=p16$Yutm)
utms = SpatialPoints(coord_utm, proj4string=CRS('+proj=utm +zone=22'))
longlats <- spTransform(utms, CRS('+init=epsg:4326'))
longlats$alt = p16$altitude
mapview(longlats, zcol = c("alt"), legend = TRUE)

Ajout de l’altitude au modèle (régression sur les données de Paracou)

(Attention ici : passage du carré de 500x500 au carré de 5x5)

Nous souhaitons rajouter l’altitude comme variable du modèle ; non pas comme un champ aléatoire mais comme un champ fixé. La syntaxe est donc légèrement différente, la fonction “f” (qu’on nommera ici “alti” comme altitude) prendra en argument model = “linear”.

Toujours dans cette même fonction, l’argument map est cette fois défini comme une seconde fonction qui retourne l’altitude à chaque coordonnée des points considérés. Cette seconde fonction nécessite un SpatialPixelsDataFrame de l’altitude, c’est à dire des valeurs d’altitude placées sur une grille régulière, ce qui n’est actuellement pas le cas puisque nos données correspondent à l’emplacement de chaque arbre. Nous devons donc effectuer au préalable un krigeage, détaillé ci-dessous.

Attention : si vous obtenez une erreur de type “Error in match.names(clabs, names(xi)) : names do not match previous names”, essayez de changer le nom de la fonction “f”… tous les noms ne sont pas acceptés et les messages d’erreur ne sont pas très clairs.

# Préparation des données
altitude = SpatialPointsDataFrame(coords = data.frame(x=p16$X/100,y=p16$Y/100), 
                                  data = data.frame(z=p16$altitude), proj4string = CRS(as.character(NA)))
sample_altitude =  SpatialPointsDataFrame(coords = data.frame(x=sample_p16$X/100,y=sample_p16$Y/100), 
                                          data = data.frame(z=sample_p16$altitude), proj4string = CRS(as.character(NA)))

sample_altitude <- sample_altitude[-zerodist(sample_altitude)[,1],]          # enlève les points de même coordonnées  


# Krigeage
vgmEmpirique <- gstat::variogram(z ~ 1, data = sample_altitude)    
vgmX <- fit.variogram(vgmEmpirique, vgm("Gau"))                              # Ajustement d'un modèle gaussien 
geoX <- gstat(formula = z ~ 1, locations = sample_altitude, model = vgmX)    # Objet geostat qui décrit toutes les                                                                                               # caractéristiques de la modélisation.

grid <- expand.grid((0:100)/20, (0:100)/20)         # Krigeage sur une grille de résolution 5cm (carré de 5m x 5m).
names(grid) <- c("x", "y")
gridded(grid) <- ~x + y
altitude_grid <- predict(geoX, newdata = grid)      # Calcul de la valeur de l'altitude sur les points de la grille       

## [using ordinary kriging]

altitude_grid=altitude_grid[,-2]
names(altitude_grid) = c("z")


# Plot
image(altitude_grid, col = topo.colors(20, alpha = 1), asp = 1)
contour(altitude_grid, add = TRUE)
title(main = "Krigeage de l'altitude", font.main = 4)

Seconde fonction qui retourne l’altitude des points étudiés

fct_alt <- function(x,y) {                              # /!\ altitude_grid DOIT ABSOLUMENT être un SpatialPixelsDataFrame
  
  spp <- SpatialPoints(data.frame(x=x,y=y))             # SpatialPoint object
  proj4string(spp) <- CRS(proj4string(altitude_grid))   # auquel on attache le système de coordonnées de l'objet altitude_grid
  v <- over(spp, altitude_grid)                         # extrait l'altitude correspondante aux coordonnées des points
  v[is.na(v)] <- 0                                      # enlève les NA
  return(v$z)
} 

Nous testons ensuite notre modèle comprenant l’altitude sur quelques espèces d’arbres.

# Choix de l'espèce et du genre 
genre = "Vouacapoua"
espece = "americana"

# Préparation des données
voam = preparation(data.frame(x = p16$X[p16$Genus==genre&p16$Species==espece]/100, y = p16$Y[p16$Genus==genre&p16$Species==espece]/100))


# Altitude seule 
fit_alt_voam = lgcp(coordinates ~ alti(map = fct_alt(x,y), model = "linear") + Intercept, data = voam$coord_pts, samplers = bnd)

lambda_alt_voam <- predict(fit_alt_voam, pix, ~ exp(alti + Intercept))

# Champ SPDE seul
model_matern = inla.spde2.pcmatern(voam$mesh, prior.sigma = c(0.1, 0.01), prior.range = c(5, 0.01))
fit_spde_voam = lgcp(coordinates ~ field(map = coordinates, model = model_matern) + Intercept, data = voam$coord_pts, samplers = bnd)

lambda_spde_voam <- predict(fit_spde_voam, pix, ~ exp(field + Intercept))

# Altitude + champ SPDE
model_matern = inla.spde2.pcmatern(voam$mesh, prior.sigma = c(0.1, 0.01), prior.range = c(5, 0.01))
fit_alt_spde_voam = lgcp(coordinates ~  field(map = coordinates, model = model_matern) + alti(map = fct_alt(x,y), model = "linear") + Intercept, data = voam$coord_pts, samplers = bnd)

lambda_alt_spde_voam <- predict(fit_alt_spde_voam, pix, ~ exp(field + alti + Intercept))
lambda_alt_spde_voam_df = as.data.frame(lambda_alt_spde_voam)


# Plot

# Semis de points
ggplot() + gg(voam$coord_pts) + gg(bnd) + coord_fixed()

## Regions defined for each Polygons

# Points avec altitude
abc = data.frame(a=p16$X[p16$Genus==genre&p16$Species==espece]/100,b=p16$Y[p16$Genus==genre&p16$Species==espece]/100,c=p16$altitude[p16$Genus==genre&p16$Species==espece])
ggplot(abc, aes(x=a, y=b)) + geom_point(aes(color = c)) + scale_color_viridis(option="D") + coord_fixed() + ggtitle(paste(genre, espece))

# Intensité moyenne prédite des 3 modèles
max_voam = max(max(lambda_alt_voam$mean),max(lambda_spde_voam$mean),max(lambda_alt_spde_voam$mean))  # Max. échelle de couleurs
ggplot() + gg(lambda_alt_voam) + gg(bnd) + gg(voam$coord_pts, shape="+") + coord_equal() + ggtitle("Altitude") + colsc(max_voam)

## Regions defined for each Polygons

## Scale for 'fill' is already present. Adding another scale for 'fill',
## which will replace the existing scale.

ggplot() + gg(lambda_spde_voam) + gg(bnd) + gg(voam$coord_pts, shape="+") + coord_equal() + ggtitle("SPDE") + colsc(max_voam)

## Regions defined for each Polygons
## Scale for 'fill' is already present. Adding another scale for 'fill',
## which will replace the existing scale.

ggplot() + gg(lambda_alt_spde_voam) + gg(bnd) + gg(voam$coord_pts, shape="+") + coord_equal() + ggtitle("Altitude + SPDE") + colsc(max_voam)

## Regions defined for each Polygons
## Scale for 'fill' is already present. Adding another scale for 'fill',
## which will replace the existing scale.

# Loi à posteriori des coefficients linéaires de l'altitude
plot(fit_alt_voam, "alti")

plot(fit_alt_spde_voam, "alti")

# Critères DIC et WAIC 
scores <- data.frame(
  dic=c(alt = fit_alt_voam$waic$waic, spde = fit_spde_voam$waic$waic, alt_spde = fit_alt_spde_voam$waic$waic),
  waic=c(alt = fit_alt_voam$dic$dic, spde = fit_spde_voam$dic$dic, alt_spde = fit_alt_spde_voam$dic$dic))
rownames(scores) <- c("altitude", "field", "altitude + field")
scores

##                        dic      waic
## altitude         -512.9292 -515.3266
## field            -566.1474 -566.5644
## altitude + field -570.1489 -570.5601

summary(fit_alt_voam)

## 
## --- Likelihoods ----------------------------------------------------------------------------------
## 
## 
## --- Criteria -------------------------------------------------------------------------------------
## 
## Watanabe-Akaike information criterion (WAIC):    -5.129e+02
## Deviance Information Criterion (DIC):        -5.153e+02
## 
## --- Fixed effects -------------------------------------------------------------------------------- 
## 
##                 mean         sd 0.025quant   0.5quant  0.975quant
## alti       0.1982447 0.01987555  0.1597474  0.1980681 0.237728367
## Intercept -0.6415450 0.33673873 -1.3205301 -0.6352843 0.002163823
##                 mode signif
## alti       0.1977187   TRUE
## Intercept -0.6228578  FALSE
## 
## --- Random effects ------------------------------------------------------------------------------- 
## 
## None.
## 
## --- All hyper parameters (internal representation) ----------------------------------------------- 
## 
## NULL

summary(fit_spde_voam)

## 
## --- Likelihoods ----------------------------------------------------------------------------------
## 
## 
## --- Criteria -------------------------------------------------------------------------------------
## 
## Watanabe-Akaike information criterion (WAIC):    -5.661e+02
## Deviance Information Criterion (DIC):        -5.666e+02
## 
## --- Fixed effects -------------------------------------------------------------------------------- 
## 
##                mean        sd 0.025quant  0.5quant 0.975quant      mode
## Intercept 0.8946979 0.6157455 -0.3557254 0.9033231   2.099715 0.9210275
##           signif
## Intercept  FALSE
## 
## --- Random effects ------------------------------------------------------------------------------- 
## 
## field ranges: mean = [-1.82194, 3.13234], sd = [0.641126, 1.05395], quantiles = [-3.96713 : 4.55821]
## 
## --- All hyper parameters (internal representation) ----------------------------------------------- 
## 
##                      mean         sd 0.025quant  0.5quant 0.975quant
## Range for field 4.4804277 0.91161583  2.9475646 4.3932478  6.5215763
## Stdev for field 0.5897562 0.08069693  0.4479424 0.5840389  0.7643381
##                      mode
## Range for field 4.2259591
## Stdev for field 0.5726136
## 
## 
## --- Field 'field' transformed hyper parameters ---
##          param       mean         var        sd         lq     median
## 1        range  4.4792521 0.818602944 0.9047668  2.9552339  4.3917884
## 2    log.range  1.4793818 0.040677686 0.2016871  1.0819332  1.4794914
## 3     variance  0.3541593 0.009519242 0.0975666  0.2010395  0.3409422
## 4 log.variance -1.0746004 0.074010033 0.2720479 -1.6062641 -1.0763325
##           uq
## 1  6.4953877
## 2  1.8723103
## 3  0.5817151
## 4 -0.5400997

summary(fit_alt_spde_voam)

## 
## --- Likelihoods ----------------------------------------------------------------------------------
## 
## 
## --- Criteria -------------------------------------------------------------------------------------
## 
## Watanabe-Akaike information criterion (WAIC):    -5.701e+02
## Deviance Information Criterion (DIC):        -5.706e+02
## 
## --- Fixed effects -------------------------------------------------------------------------------- 
## 
##                 mean         sd  0.025quant   0.5quant 0.975quant
## alti       0.1152663 0.02337132  0.06936644  0.1152541  0.1611970
## Intercept -0.5149656 0.67676818 -1.86962484 -0.5108825  0.8168936
##                 mode signif
## alti       0.1152328   TRUE
## Intercept -0.5020177  FALSE
## 
## --- Random effects ------------------------------------------------------------------------------- 
## 
## field ranges: mean = [-1.7757, 2.37553], sd = [0.620763, 0.894035], quantiles = [-3.61457 : 3.76245]
## 
## --- All hyper parameters (internal representation) ----------------------------------------------- 
## 
##                      mean       sd 0.025quant  0.5quant 0.975quant
## Range for field 5.2571960 1.168258   3.330659 5.1333402  7.9050183
## Stdev for field 0.4874612 0.075589   0.355482 0.4817944  0.6522584
##                      mode
## Range for field 4.8961122
## Stdev for field 0.4708139
## 
## 
## --- Field 'field' transformed hyper parameters ---
##          param       mean         var         sd         lq     median
## 1        range  5.2555776 1.343533145 1.15910877  3.3382363  5.1314740
## 2    log.range  1.6354701 0.048247758 0.21965372  1.2036900  1.6351308
## 3     variance  0.2431839 0.005784601 0.07605656  0.1268349  0.2319968
## 4 log.variance -1.4608678 0.095105293 0.30839146 -2.0670681 -1.4613436
##           uq
## 1  7.8707730
## 2  2.0644874
## 3  0.4230275
## 4 -0.8584636

# Intensité observée avec density
densi = density(voam$coord_pts_ppp, sigma = bw.diggle(voam$coord_pts_ppp))
densi$v[densi$v<0]=0

# Intensité moyenne prédite
mat_int = matrix(lambda_alt_spde_voam_df$mean, ncol=200, nrow=200, byrow = TRUE)

# Ecart-type de la prédiction
mat_eca = matrix(lambda_alt_spde_voam_df$sd, ncol=200, nrow=200, byrow = TRUE)

# Echelle de couleurs commune pour intensité obs. et moy. prédite
intensity_max <- ceiling(max(max(densi), max(lambda_alt_spde_voam_df$mean)))
col_fixe <- colourmap(hcl.colors(intensity_max), breaks=0:intensity_max)

# Plot
plot(densi, main="Intensité observée", col=col_fixe)

plot(im(mat_int, xcol=seq(ncol(mat_int)), yrow=seq(nrow(mat_int))), main="Intensité moy. prédite", col=col_fixe)

plot(im(mat_eca, xcol=seq(ncol(mat_eca)), yrow=seq(nrow(mat_eca))), main="Ecart-type") 

post.range = spde.posterior(fit_alt_spde_voam, name="field", what="range"); plot(post.range)

post.matcorr = spde.posterior(fit_alt_spde_voam, name="field",what="matern.correlation"); plot(post.matcorr)

# Choix de l'espèce et du genre 
genre = "Eperua"
espece = "falcata"

# Préparation des données
epfa = preparation(data.frame(x = p16$X[p16$Genus==genre&p16$Species==espece]/100, y = p16$Y[p16$Genus==genre&p16$Species==espece]/100))


# Altitude seule 
fit_alt_epfa = lgcp(coordinates ~ alti(map = fct_alt(x,y), model = "linear") + Intercept, data = epfa$coord_pts, samplers = bnd)

lambda_alt_epfa <- predict(fit_alt_epfa, pix, ~ exp(alti + Intercept))

# Champ SPDE seul
model_matern = inla.spde2.pcmatern(epfa$mesh, prior.sigma = c(0.1, 0.01), prior.range = c(5, 0.01))
fit_spde_epfa = lgcp(coordinates ~ field(map = coordinates, model = model_matern) + Intercept, data = epfa$coord_pts, samplers = bnd)

lambda_spde_epfa <- predict(fit_spde_epfa, pix, ~ exp(field + Intercept))

# Altitude + champ SPDE
model_matern = inla.spde2.pcmatern(epfa$mesh, prior.sigma = c(0.1, 0.01), prior.range = c(5, 0.01))
fit_alt_spde_epfa = lgcp(coordinates ~  field(map = coordinates, model = model_matern) + alti(map = fct_alt(x,y), model = "linear") + Intercept, data = epfa$coord_pts, samplers = bnd)

lambda_alt_spde_epfa <- predict(fit_alt_spde_epfa, pix, ~ exp(field + alti + Intercept))
lambda_alt_spde_df_epfa = as.data.frame(lambda_alt_spde_epfa)


# Plot

# Semis de points
ggplot() + gg(epfa$coord_pts) + gg(bnd) + coord_fixed()

## Regions defined for each Polygons

# Points avec altitude
abc = data.frame(a=p16$X[p16$Genus==genre&p16$Species==espece]/100,b=p16$Y[p16$Genus==genre&p16$Species==espece]/100,c=p16$altitude[p16$Genus==genre&p16$Species==espece])
ggplot(abc, aes(x=a, y=b)) + geom_point(aes(color = c)) + scale_color_viridis(option="D") + coord_fixed() + ggtitle(paste(genre, espece))

# Intensité moyenne prédite des 3 modèles
max_epfa = max(max(lambda_alt_epfa$mean),max(lambda_spde_epfa$mean),max(lambda_alt_spde_epfa$mean))  # max. échelle de couleurs
ggplot() + gg(lambda_alt_epfa) + gg(bnd) + gg(epfa$coord_pts, shape="+") + coord_equal() + ggtitle("Altitude") + colsc(max_epfa)

## Regions defined for each Polygons

## Scale for 'fill' is already present. Adding another scale for 'fill',
## which will replace the existing scale.

ggplot() + gg(lambda_spde_epfa) + gg(bnd) + gg(epfa$coord_pts, shape="+") + coord_equal() + ggtitle("SPDE") + colsc(max_epfa)

## Regions defined for each Polygons
## Scale for 'fill' is already present. Adding another scale for 'fill',
## which will replace the existing scale.

ggplot() + gg(lambda_alt_spde_epfa) + gg(bnd) + gg(epfa$coord_pts, shape="+") + coord_equal() + ggtitle("Altitude + SPDE") + colsc(max_epfa)

## Regions defined for each Polygons
## Scale for 'fill' is already present. Adding another scale for 'fill',
## which will replace the existing scale.

# Loi à posteriori des coefficients linéaires de l'altitude
plot(fit_alt_epfa, "alti")

plot(fit_alt_spde_epfa, "alti")

# Critères DIC et WAIC
scores <- data.frame(
  dic=c(alt = fit_alt_epfa$waic$waic, spde = fit_spde_epfa$waic$waic, alt_spde = fit_alt_spde_epfa$waic$waic),
  waic=c(alt = fit_alt_epfa$dic$dic, spde = fit_spde_epfa$dic$dic, alt_spde = fit_alt_spde_epfa$dic$dic))
rownames(scores) <- c("altitude", "field", "altitude + field")
scores

##                        dic      waic
## altitude         -143.1317 -143.7463
## field            -609.1629 -607.6059
## altitude + field -608.5919 -606.8954

summary(fit_alt_epfa)

## 
## --- Likelihoods ----------------------------------------------------------------------------------
## 
## 
## --- Criteria -------------------------------------------------------------------------------------
## 
## Watanabe-Akaike information criterion (WAIC):    -1.431e+02
## Deviance Information Criterion (DIC):        -1.437e+02
## 
## --- Fixed effects -------------------------------------------------------------------------------- 
## 
##                  mean         sd 0.025quant   0.5quant  0.975quant
## alti      -0.06904678 0.02043147 -0.1091926 -0.0690363 -0.02899805
## Intercept  2.41347156 0.23360976  1.9413661  2.4182071  2.85894137
##                  mode signif
## alti      -0.06901356   TRUE
## Intercept  2.42760840   TRUE
## 
## --- Random effects ------------------------------------------------------------------------------- 
## 
## None.
## 
## --- All hyper parameters (internal representation) ----------------------------------------------- 
## 
## NULL

summary(fit_spde_epfa)

## 
## --- Likelihoods ----------------------------------------------------------------------------------
## 
## 
## --- Criteria -------------------------------------------------------------------------------------
## 
## Watanabe-Akaike information criterion (WAIC):    -6.092e+02
## Deviance Information Criterion (DIC):        -6.076e+02
## 
## --- Fixed effects -------------------------------------------------------------------------------- 
## 
##                 mean        sd 0.025quant   0.5quant 0.975quant       mode
## Intercept -0.3502781 0.6574403  -1.697013 -0.3326241   0.899974 -0.2975257
##           signif
## Intercept  FALSE
## 
## --- Random effects ------------------------------------------------------------------------------- 
## 
## field ranges: mean = [-1.77624, 5.74104], sd = [0.717743, 2.02388], quantiles = [-5.4864 : 8.35587]
## 
## --- All hyper parameters (internal representation) ----------------------------------------------- 
## 
##                     mean        sd 0.025quant 0.5quant 0.975quant     mode
## Range for field 2.315757 0.3129724  1.7667776 2.293108   2.992918 2.247193
## Stdev for field 1.142503 0.1078816  0.9458953 1.137194   1.369543 1.126442
## 
## 
## --- Field 'field' transformed hyper parameters ---
##          param      mean        var        sd         lq    median
## 1        range 2.3155014 0.09664505 0.3108779  1.7686429 2.2926387
## 2    log.range 0.8307180 0.01802728 0.1342657  0.5690927 0.8295366
## 3     variance 1.3166633 0.06183539 0.2486672  0.8963457 1.2928225
## 4 log.variance 0.2575794 0.03545913 0.1883059 -0.1109539 0.2566033
##          uq
## 1 2.9874274
## 2 1.0952718
## 3 1.8698500
## 4 0.6270328

summary(fit_alt_spde_epfa)

## 
## --- Likelihoods ----------------------------------------------------------------------------------
## 
## 
## --- Criteria -------------------------------------------------------------------------------------
## 
## Watanabe-Akaike information criterion (WAIC):    -6.086e+02
## Deviance Information Criterion (DIC):        -6.069e+02
## 
## --- Fixed effects -------------------------------------------------------------------------------- 
## 
##                  mean         sd 0.025quant    0.5quant 0.975quant
## alti      -0.03269948 0.05442809 -0.1408992 -0.03231621 0.07325973
## Intercept  0.03650594 0.91217420 -1.8273892  0.06175723 1.75923485
##                  mode signif
## alti      -0.03156701  FALSE
## Intercept  0.11160607  FALSE
## 
## --- Random effects ------------------------------------------------------------------------------- 
## 
## field ranges: mean = [-1.73656, 5.70725], sd = [0.721576, 2.03381], quantiles = [-5.43981 : 8.25936]
## 
## --- All hyper parameters (internal representation) ----------------------------------------------- 
## 
##                     mean        sd 0.025quant 0.5quant 0.975quant     mode
## Range for field 2.274916 0.3138912  1.7250694 2.252120   2.954831 2.206272
## Stdev for field 1.143510 0.1087392  0.9440265 1.138576   1.371684 1.129071
## 
## 
## --- Field 'field' transformed hyper parameters ---
##          param     mean        var        sd         lq    median
## 1        range 2.274655 0.09721038 0.3117858  1.7266341 2.2516470
## 2    log.range 0.812541 0.01879556 0.1370969  0.5450284 0.8114918
## 3     variance 1.319140 0.06284580 0.2506907  0.8934717 1.2959568
## 4 log.variance 0.259194 0.03602727 0.1898085 -0.1141936 0.2590217
##          uq
## 1 2.9489123
## 2 1.0823099
## 3 1.8748176
## 4 0.6296763

# Intensité observée avec density
densi = density(epfa$coord_pts_ppp, sigma = bw.diggle(epfa$coord_pts_ppp)) 
densi$v[densi$v<0]=0


# Intensité moyenne prédite
mat_int = matrix(lambda_alt_spde_df_epfa$mean, ncol=200, nrow=200, byrow = TRUE)

# Ecart-type de la prédiction
mat_eca = matrix(lambda_alt_spde_df_epfa$sd, ncol=200, nrow=200, byrow = TRUE)    


# Echelle de couleurs commune pour l'intensité obs. et moy. prédite
intensity_max <- ceiling(max(max(densi), max(lambda_alt_spde_df_epfa$mean)))
col_fixe <- colourmap(hcl.colors(intensity_max), breaks=0:intensity_max)

# Plot
plot(densi, main="Intensité observée", col=col_fixe)

plot(im(mat_int, xcol=seq(ncol(mat_int)), yrow=seq(nrow(mat_int))), main="Intensité moy. prédite", col=col_fixe)

plot(im(mat_eca, xcol=seq(ncol(mat_eca)), yrow=seq(nrow(mat_eca))), main="Ecart-type") 

post.range = spde.posterior(fit_alt_spde_epfa, name="field", what="range"); plot(post.range)

post.matcorr = spde.posterior(fit_alt_spde_epfa, name="field",what="matern.correlation"); plot(post.matcorr)

# Choix de l'espèce et du genre 
genre = "Oenocarpus" 
espece = "bataua"

# Préparation des données
oeba = preparation(data.frame(x = p16$X[p16$Genus==genre&p16$Species==espece]/100, y = p16$Y[p16$Genus==genre&p16$Species==espece]/100))


# Altitude seule 
fit_alt_oeba = lgcp(coordinates ~ alti(map = fct_alt(x,y), model = "linear") + Intercept, data = oeba$coord_pts, samplers = bnd)

lambda_alt_oeba <- predict(fit_alt_oeba, pix, ~ exp(alti + Intercept))

# Champ SPDE seul
model_matern = inla.spde2.pcmatern(oeba$mesh, prior.sigma = c(0.1, 0.01), prior.range = c(5, 0.01))
fit_spde_oeba = lgcp(coordinates ~ field(map = coordinates, model = model_matern) + Intercept, data = oeba$coord_pts, samplers = bnd)

lambda_spde_oeba <- predict(fit_spde_oeba, pix, ~ exp(field + Intercept))

# Altitude + champ SPDE
model_matern = inla.spde2.pcmatern(oeba$mesh, prior.sigma = c(0.1, 0.01), prior.range = c(5, 0.01))
fit_alt_spde_oeba = lgcp(coordinates ~  field(map = coordinates, model = model_matern) + alti(map = fct_alt(x,y), model = "linear") + Intercept, data = oeba$coord_pts, samplers = bnd)

lambda_alt_spde_oeba <- predict(fit_alt_spde_oeba, pix, ~ exp(field + alti + Intercept))
lambda_alt_spde_df_oeba = as.data.frame(lambda_alt_spde_oeba)


# Plot

# Semis de points
ggplot() + gg(oeba$coord_pts) + gg(bnd) + coord_fixed()

## Regions defined for each Polygons

# Points avec altitude
abc = data.frame(a=p16$X[p16$Genus==genre&p16$Species==espece]/100,b=p16$Y[p16$Genus==genre&p16$Species==espece]/100,c=p16$altitude[p16$Genus==genre&p16$Species==espece])
ggplot(abc, aes(x=a, y=b)) + geom_point(aes(color = c)) + scale_color_viridis(option="D") + coord_fixed() + ggtitle(paste(genre, espece))

# Intensité moyenne prédite des 3 modèles
max_oeba = max(max(lambda_alt_oeba$mean),max(lambda_spde_oeba$mean),max(lambda_alt_spde_oeba$mean)) # max. échelle de couleurs
ggplot() + gg(lambda_alt_oeba) + gg(bnd) + gg(oeba$coord_pts, shape="+") + coord_equal() + ggtitle("Altitude") + colsc(max_oeba)

## Regions defined for each Polygons

## Scale for 'fill' is already present. Adding another scale for 'fill',
## which will replace the existing scale.

ggplot() + gg(lambda_spde_oeba) + gg(bnd) + gg(oeba$coord_pts, shape="+") + coord_equal() + ggtitle("SPDE") + colsc(max_oeba)

## Regions defined for each Polygons
## Scale for 'fill' is already present. Adding another scale for 'fill',
## which will replace the existing scale.

ggplot() + gg(lambda_alt_spde_oeba) + gg(bnd) + gg(oeba$coord_pts, shape="+") + coord_equal() + ggtitle("Altitude + SPDE") + colsc(max_oeba)

## Regions defined for each Polygons
## Scale for 'fill' is already present. Adding another scale for 'fill',
## which will replace the existing scale.

# Loi à posteriori des coefficients linéaires de l'altitude
plot(fit_alt_oeba, "alti")

plot(fit_alt_spde_oeba, "alti")

# Critères DIC et WAIC
scores <- data.frame(
  dic=c(alt = fit_alt_oeba$waic$waic, spde = fit_spde_oeba$waic$waic, alt_spde = fit_alt_spde_oeba$waic$waic),
  waic=c(alt = fit_alt_oeba$dic$dic, spde = fit_spde_oeba$dic$dic, alt_spde = fit_alt_spde_oeba$dic$dic))
rownames(scores) <- c("altitude", "field", "altitude + field")
scores

##                        dic      waic
## altitude         -901.5785 -903.5977
## field            -933.0700 -933.6494
## altitude + field -933.2780 -933.7714

summary(fit_alt_oeba)

## 
## --- Likelihoods ----------------------------------------------------------------------------------
## 
## 
## --- Criteria -------------------------------------------------------------------------------------
## 
## Watanabe-Akaike information criterion (WAIC):    -9.016e+02
## Deviance Information Criterion (DIC):        -9.036e+02
## 
## --- Fixed effects -------------------------------------------------------------------------------- 
## 
##                  mean         sd  0.025quant    0.5quant 0.975quant
## alti      -0.01179511 0.01206122 -0.03544741 -0.01180521 0.01189162
## Intercept  2.64143243 0.15674084  2.32830073  2.64329762 2.94393736
##                  mode signif
## alti      -0.01182442  FALSE
## Intercept  2.64700462   TRUE
## 
## --- Random effects ------------------------------------------------------------------------------- 
## 
## None.
## 
## --- All hyper parameters (internal representation) ----------------------------------------------- 
## 
## NULL

summary(fit_spde_oeba)

## 
## --- Likelihoods ----------------------------------------------------------------------------------
## 
## 
## --- Criteria -------------------------------------------------------------------------------------
## 
## Watanabe-Akaike information criterion (WAIC):    -9.331e+02
## Deviance Information Criterion (DIC):        -9.336e+02
## 
## --- Fixed effects -------------------------------------------------------------------------------- 
## 
##               mean        sd 0.025quant 0.5quant 0.975quant    mode signif
## Intercept 2.459096 0.3387078   1.759711 2.461535   3.150623 2.46739   TRUE
## 
## --- Random effects ------------------------------------------------------------------------------- 
## 
## field ranges: mean = [-0.528934, 0.491626], sd = [0.3415, 0.411035], quantiles = [-1.40309 : 1.30928]
## 
## --- All hyper parameters (internal representation) ----------------------------------------------- 
## 
##                      mean        sd 0.025quant  0.5quant 0.975quant
## Range for field 7.9403136 2.6526225 4.08837086 7.4846036 14.3944840
## Stdev for field 0.1767393 0.0477657 0.09970707 0.1712971  0.2861291
##                      mode
## Range for field 6.6663295
## Stdev for field 0.1609388
## 
## 
## --- Field 'field' transformed hyper parameters ---
##          param        mean          var         sd           lq
## 1        range  7.93466725 6.8838889509 2.62371663  4.106581538
## 2    log.range  2.02025022 0.1022367481 0.31974482  1.410394843
## 3     variance  0.03345191 0.0003457436 0.01859418  0.009957444
## 4 log.variance -3.53770321 0.2877557107 0.53642866 -4.611511038
##        median          uq
## 1  7.47955216 14.30889014
## 2  2.01173690  2.66281840
## 3  0.02930361  0.08105774
## 4 -3.53050081 -2.50984618

summary(fit_alt_spde_oeba)

## 
## --- Likelihoods ----------------------------------------------------------------------------------
## 
## 
## --- Criteria -------------------------------------------------------------------------------------
## 
## Watanabe-Akaike information criterion (WAIC):    -9.333e+02
## Deviance Information Criterion (DIC):        -9.338e+02
## 
## --- Fixed effects -------------------------------------------------------------------------------- 
## 
##                  mean         sd  0.025quant    0.5quant  0.975quant
## alti      -0.02434751 0.01451099 -0.05299715 -0.02430024 0.004005865
## Intercept  2.76132361 0.37867946  1.99356014  2.76410644 3.514588434
##                  mode signif
## alti      -0.02420687  FALSE
## Intercept  2.76985799   TRUE
## 
## --- Random effects ------------------------------------------------------------------------------- 
## 
## field ranges: mean = [-0.558687, 0.423791], sd = [0.337832, 0.415136], quantiles = [-1.44352 : 1.22923]
## 
## --- All hyper parameters (internal representation) ----------------------------------------------- 
## 
##                      mean         sd 0.025quant  0.5quant 0.975quant
## Range for field 7.9407660 2.67227909 4.06742403 7.4796476 14.4480425
## Stdev for field 0.1745113 0.04814275 0.09698259 0.1689933  0.2847932
##                      mode
## Range for field 6.6529763
## Stdev for field 0.1584815
## 
## 
## --- Field 'field' transformed hyper parameters ---
##          param        mean          var         sd          lq      median
## 1        range  7.93503961 6.9853394734 2.64297928  4.08581952  7.47456187
## 2    log.range  2.01958195 0.1036717215 0.32198093  1.40531969  2.01106915
## 3     variance  0.03270457 0.0003445674 0.01856253  0.00943151  0.02851929
## 4 log.variance -3.56611774 0.3003139629 0.54800909 -4.66568056 -3.55762651
##            uq
## 1 14.36165942
## 2  2.66651263
## 3  0.08033236
## 4 -2.51880418

# Intensité observée avec density
densi = density(oeba$coord_pts_ppp, sigma = bw.diggle(oeba$coord_pts_ppp)) 
densi$v[densi$v<0]=0

# Intensité moyenne prédite
mat_int = matrix(lambda_alt_spde_df_oeba$mean, ncol=200, nrow=200, byrow = TRUE)

# Ecart-type de la prédiction
mat_eca = matrix(lambda_alt_spde_df_oeba$sd, ncol=200, nrow=200, byrow = TRUE)     



# Echelle de couleurs commune pour l'intensité obs. et moy. prédite
intensity_max <- ceiling(max(max(densi), max(lambda_alt_spde_df_oeba$mean)))
col_fixe <- colourmap(hcl.colors(intensity_max), breaks=0:intensity_max)

# Plot
plot(densi, main="Intensité observée", col=col_fixe)

plot(im(mat_int, xcol=seq(ncol(mat_int)), yrow=seq(nrow(mat_int))), main="Intensité moy. prédite", col=col_fixe)

plot(im(mat_eca, xcol=seq(ncol(mat_eca)), yrow=seq(nrow(mat_eca))), main="Ecart-type")

post.range = spde.posterior(fit_alt_spde_oeba, name="field", what="range"); plot(post.range)

post.matcorr = spde.posterior(fit_alt_spde_oeba, name="field",what="matern.correlation"); plot(post.matcorr)

Qualea rosea et Vouacapoua americana : deux espèces qui se repoussent.

# Préparation des deux espèces

genre = "Qualea"
espece = "rosea"
quro_grd = preparation(data.frame(x = p16$X[p16$Genus==genre&p16$Species==espece], y = p16$Y[p16$Genus==genre&p16$Species==espece]), size=500)

genre = "Vouacapoua"
espece = "americana"
voam_grd = preparation(data.frame(x = p16$X[p16$Genus==genre&p16$Species==espece], y = p16$Y[p16$Genus==genre&p16$Species==espece]), size=500)


# Graphique des deux espèces

doubleplot = data.frame( rbind(quro_grd$df,voam_grd$df) , Espece = c(rep("Qualea rosea",quro_grd$nb_obs),rep("Vouacapoua americana",voam_grd$nb_obs)) )
ggplot(doubleplot, aes(x=x, y=y)) + geom_point(aes(color = Espece)) + coord_fixed() 

# Calcul de la fonction M

X = c(voam_grd$df[,1],quro_grd$df[,1])
Y = c(voam_grd$df[,2],quro_grd$df[,2])
PointType <- c(rep("V. Americana", voam_grd$nb_obs), rep("Q. Rosea", quro_grd$nb_obs))
PointWeight <- rep(1,voam_grd$nb_obs+quro_grd$nb_obs)
pattern <- wmppp(data.frame(X, Y, PointType, PointWeight), owin(xrange=c(0,500), yrange=c(0,500)))


result_M = Mhat(X=pattern, ReferenceType="V. Americana", NeighborType="Q. Rosea")
plot(result_M$r, result_M$M, type = "l", ylim=c(0.15,1), xlab="r", ylab="M(r)")
lines(result_M$r, result_M$theo, type = "l", col = "red", lty=2)
legend(x=250, y=0.32, legend=c("M indépendant", "M observé"), col=c("red", "black"), lty=c(2,1))

Cas homogène

# Processus n°1 : Poisson homogène (lambda=20) sur la moitié du domaine
poish_moitie = rpoispp(lambda = 20, win = owin(xrange=c(0,2.5), yrange=c(0,5)))
processus_1 = preparation(data.frame(x=poish_moitie$x, y= poish_moitie$y))

# Processus n°2 : processus n°1 auquel on a rajouté deux points aux deux coins opposés du carré (de coord. (5,0) et (5,5)). 
poish_moitie_deux_coins = data.frame(x = c(processus_1$df$x, 5, 5), y= c(processus_1$df$y, 0, 5))
processus_2 = preparation(data.frame(x=poish_moitie_deux_coins$x, y= poish_moitie_deux_coins$y))

# Processus n°3 : processus n°2 mais dont les points ne sont cette fois plus sur la frontière mais légèrement rentrés dans le domaine (de coord. (4.8,0.2) et (4.8,4.8)). 
poish_moitie_deux_coins_bis = data.frame(x = c(processus_1$df$x, 4.8, 4.8), y= c(processus_1$df$y, 0.2, 4.8))
processus_3 = preparation(data.frame(x=poish_moitie_deux_coins_bis$x, y= poish_moitie_deux_coins_bis$y))


# Régression de Cox homogène
fit_coxh1 <- lgcp(coordinates ~ Intercept, data = processus_1$coord_pts, samplers = bnd)

fit_coxh2 <- lgcp(coordinates ~ Intercept, data = processus_2$coord_pts, samplers = bnd)

fit_coxh3 <- lgcp(coordinates ~ Intercept, data = processus_2$coord_pts, samplers = bnd)

# Plot
ggplot() + gg(processus_1$coord_pts) + gg(bnd) + coord_fixed()

## Regions defined for each Polygons

ggplot() + gg(processus_2$coord_pts) + gg(bnd) + coord_fixed()

## Regions defined for each Polygons

ggplot() + gg(processus_3$coord_pts) + gg(bnd) + coord_fixed()

## Regions defined for each Polygons

# Intensités
print(paste("Intensité observée processus n°1 = ", (processus_1$nb_obs)/25))             

## [1] "Intensité observée processus n°1 =  10.88"

print(paste("Intensité cox homogène processus n°1 = ", exp(fit_coxh1$summary.fixed$mean)))

## [1] "Intensité cox homogène processus n°1 =  19.8342735787074"

print(paste("Intensité observée processus n°2 = ", (processus_2$nb_obs)/25))             

## [1] "Intensité observée processus n°2 =  10.96"

print(paste("Intensité cox homogène processus n°2 = ", exp(fit_coxh2$summary.fixed$mean)))

## [1] "Intensité cox homogène processus n°2 =  10.9596934073662"

print(paste("Intensité observée processus n°3 = ", (processus_3$nb_obs+2)/25))             

## [1] "Intensité observée processus n°3 =  11.04"

print(paste("Intensité cox homogène processus n°3 = ", exp(fit_coxh3$summary.fixed$mean)))

## [1] "Intensité cox homogène processus n°3 =  10.9596934073662"

L’intensité cox homogène processus n°1 est généralement légèrement inférieure à 20 (intensité du poisson sur la moitié gauche du domaine) : la fonction ne semble pas prendre en considération les zones vides du domaine.

Cas non-homogène

# Régression de Cox non-homogène processus n°1
model_matern = inla.spde2.pcmatern(processus_1$mesh, prior.sigma = c(0.1, 0.01), prior.range = c(5, 0.01))
fit_cox1 = lgcp(coordinates ~ field(map = coordinates, model = model_matern) + Intercept, data = processus_1$coord_pts, samplers = bnd)

lambda_cox1 <- predict(fit_cox1, pix, ~ exp(field + Intercept))

# Régression de Cox non-homogène processus n°2
model_matern = inla.spde2.pcmatern(processus_2$mesh, prior.sigma = c(0.1, 0.01), prior.range = c(5, 0.01))
fit_cox2 = lgcp(coordinates ~ field(map = coordinates, model = model_matern) + Intercept, data = processus_2$coord_pts, samplers = bnd)

lambda_cox2 <- predict(fit_cox2, pix, ~ exp(field + Intercept))

# Régression de Cox non-homogène processus n°3
model_matern = inla.spde2.pcmatern(processus_3$mesh, prior.sigma = c(0.1, 0.01), prior.range = c(5, 0.01))
fit_cox3 = lgcp(coordinates ~ field(map = coordinates, model = model_matern) + Intercept, data = processus_3$coord_pts, samplers = bnd)

lambda_cox3 <- predict(fit_cox3, pix, ~ exp(field + Intercept))


# Plot
plot(lambda_cox1)

plot(lambda_cox2)

plot(lambda_cox3)

Le problème d’implémentation n’impacte cependant pas la régression lorsque celle-ci est non-homogène.